arXiv:1502.02432v3 [physics.plasm-ph] 18 Apr 2016 


Classical Radiation Reaction in Particle-In-Cell Simulations 


M. Vranic^, J. L. Martins^, R. A. Fonseca^*’, L. O. Silva'* 

^GoLP/Instituto de Plasmas e Fusdo Nuclear, Institute Superior Tecnico, Universidade de Lisboa, Lisbon, Portugal 
^DCTI/ISCTE - Institute Universitdrio de Lisboa, 1649-026 Lisboa, Portugal 


Abstract 

Under the presence of ultra high intensity lasers or other intense electromagnetic fields the motion of particles in the 
ultrarelativistic regime can be severely affected by radiation reaction. The standard particle-in-cell (PIC) algorithms 
do not include radiation reaction effects. Even though this is a well known mechanism, there is not yet a definite 
algorithm nor a standard technique to include radiation reaction in PIC codes. We have compared several models 
for the calculation of the radiation reaction force, with the goal of implementing an algorithm for classical radiation 
reaction in the Osiris framework, a state-of-the-art PIC code. The results of the different models are compared with 
standard analytical results, and the relevance/advantages of each model are discussed. Numerical issues relevant to 
PIC codes such as resolution requirements, application of radiation reaction to macro particles and computational cost 
are also addressed. For parameters of interest where the classical description of the electron motion is applicable, all 
the models considered are shown to give comparable results. The Landau and Lifshitz reduced model is chosen for 
implementation as one of the candidates with the minimal overhead and no additional memory requirements. 
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1. Introduction 

The next generation of high-power lasers is going to 
reach intensities that will open new doors for explor¬ 
ing a wide range of physical problems with an even 
wider range of applications. The ELI project [I] is ex¬ 
pected to reach laser intensities several orders of mag¬ 
nitude higher than those available today. At intensities 
/ 10^^-10^^ W/cm^ one can expect electron-positron 

pair production mill II a. In astrophysics these inten¬ 
sities are relevant for the study of pulsars, blazars, and 
gamma-ray bursts El. High intensity laser-mater inter¬ 
actions can also produce proton and heavy ion beams 
QElilllol that are of great significance for many ap¬ 
plications, the most important being cancer treatment. 

At high intensities particle acceleration can be 
severely limited by the radiation reaction associated 
with the energy loss via radiation emission IfTTI . This 
is important whenever the radiated energy is compara¬ 
ble to the total particle energy. The threshold electro¬ 
magnetic field intensities to drive these effects vary in 
the available literature, but they are usually in the range 
10^^ - 10^^ W/cm^. Many authors have discussed the 
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classical radiation reaction in pursuit of a proper analyt¬ 
ical description lfT2l[ni[T4lfrarm[T7l[Al[T^l2nilIT1 
or experimental signatures Il22l |23] |24| |25] 123 |22l l28l 
|29l 123 laU |32l. There is also a rising interest in the 
effect of the radiation reaction on particle dynamics in 
astrophysical phenomena ll33l l34l l35l 1231 . In order to 
perform reliable particle-in-cell (PIC) simulations in the 
classical radiation-dominated regime, radiation reaction 
(RR) must be included in the equations of motion for 
the particles. Though several models for classical radia¬ 
tion reaction have been proposed in the literature, there 
is not a definite standard choice to apply in PIC codes. 
In this paper we compare several different models in or¬ 
der to find the most appropriate one for implementation 
in OSIRIS EllSSl. The chosen model should capture 
all the relevant physics in the scenarios we are aiming 
to explore at the lowest possible cost in performance. 
We test each of them with well studied examples of par¬ 
ticle motion in electromagnetic fields where the trajec¬ 
tory can be analytically expressed and estimates for the 
radiated power/energy can be obtained. Our analysis 
shows that even though there are conceptual differences 
between the models considered, in the parameters of in¬ 
terest that could be tested with near-future laser technol¬ 
ogy all the models give the same description of particle 
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motion. To observe differences, one would need to im¬ 
pose a linear acceleration so high that an extreme elec¬ 
tric field required to produce it would render a classical 
description of an electron trajectory inapplicable. The 
choice can then be made on the basis of the computa¬ 
tional overhead and with this aim the additional compu¬ 
tational cost that the implementation of each model in¬ 
troduces is presented. Among the models with the low¬ 
est computational requirements, we opted to introduce 
Landau and Lifshitz reduced model in OSIRIS frame¬ 
work. Specific questions associated with the interpre¬ 
tation of classical RR in PIC simulations are also ad¬ 
dressed. 

This paper is structured as follows. In section 2 we 
introduce the radiation reaction models and deal with 
macro particle interpretation. In section 3, the behaviour 
of all the models is investigated in the standard cases 
of synchrotron radiation and bremsstrahlung. Section 
4 underlines the difference between the results with and 
without radiation reaction for an electron in a laser pulse 
field and gives an estimate of the threshold for the de¬ 
tectable radiation reaction. The issue of the optimal 
temporal resolution is addressed in section 5, and sec¬ 
tion 6 contains estimates for the computational overhead 
for each model. Finally, in section 7 we state the con¬ 
clusions. 

2. Radiation reaction models 

The charged particle motion with radiation reaction is 
expressed by the Lorentz-Abraham-Dirac (LAD) equa¬ 
tion ||39]| : 

^ = where 

dr ^ ^ 

j,RR ^ 

^ 3mc^ \ dr^ \ dr dr JJ 

Here, denotes the electromagnetic force four- 
vector, Pfj is the particle momentum four-vector, e, m 
are the elementary charge and particle mass respec¬ 
tively, and c is the speed of light. Equation ([!]) is de¬ 
rived for a point charge. Unphysical solutions appear, 
for example, when = 0, where in addition to the 
solution with a constant velocity, eq. Q has a solution 
where the particle accelerates infinitely (the so-called 
’’runaway solution”). The principle of causality is also 
violated here with pre-acceleration solutions - these so¬ 
lutions anticipate the change of the force, so the particle 
accelerates before the force has been applied. The de¬ 
tailed explanation of these problems and suggestions for 
possible improvements are given in ll40ll4TII . 


However, even if these problems would be solved, 
the LAD equation is inconvenient for numerical inte¬ 
gration. It is possible to integrate it backwards in time 
mu, but this is not applicable to many-body problems 
ES). Therefore, the LAD equation, expressed by ([T]i, is 
inadequate for use in a simulation code. 

Numerous approximate models have been explored to 
eliminate the above-mentioned problems and to obtain 
a more efhcient computational approach. We consider 
a subset of these models appropriate for PIC implemen¬ 
tation: BOS m, LL Il44l (used also in Q and ||6l), S09 
EH, H08 EH, and F93 EH- We also consider the LL 
reduced (LLR) model BSl where the spatio-temporal 
derivatives of the fields are discarded in the equation 
of motion because they are shown to have smaller con¬ 
tribution than the particle spin, that becomes important 
only in the quantum regime. Model BOS keeps only the 
leading-order term of the LL equation. The HOS model 
estimates the total energy radiated via the Larmour for¬ 
mula, and then this energy is discounted from the parti¬ 
cle in the end of the timestep. All the models apply to 
the case of relativistic motion that is required to have an 
appreciable energy loss due to radiation emission. Their 
validity is limited to the classical domain, where the par¬ 
ticle trajectory can still be considered to be a smooth 
function of time (the individual emission events do not 
take a great fraction of the particle energy, or in other 
words, the emitted energy in the Lorentz frame momen¬ 
tarily co-moving with the emitting particle is small com¬ 
pared with mc^). Models LL and S09 also appear in Ref. 
B9l that compares the asymptotic solution for equa¬ 
tions of motion coming from strong field quantum elec¬ 
trodynamics with several classical equations of motion; 
their findings show that LL model is consistent with the 
asymptotic strong field QED description. Therefore, we 
will identify where the results obtained with other mod¬ 
els are the same (or close enough) as with LL within the 
limits of the validity of the classical description. 

To facilitate the analysis for the PIC implementation, 
we will use the 3-vector form of the equations through¬ 
out. The total change of momentum in time depends on 
the Lorentz force (F/,) and the radiation reaction force 
(F«r): 

^ = Fi + F«« (2) 

dt 

where in CGS units, is given by: 

FL = e(E + -^xB). (3) 

\ ymc I 

The radiation reaction force (Fr^) for all the considered 
models, in CGS units, is presented in Table 1. Here the 


2 





radiation back reaction is explicitly given as an addi¬ 
tional force acting on the particle, expressed as a func¬ 
tion of the electromagnetic fields E, B, the particle mo¬ 
mentum p, charge e, mass m and relativistic factor j, the 
speed of light c and time f. 

Particle-in-cell codes usually employ normalised 
units to bring all the quantities to similar orders of 
magnitude and express the physics as a function of 
fundamental plasma parameters. The normalisation in 
OSIRIS is as follows: f —> faj„, x —> xw„/c, p —> 
p/mc = yv/c, E —» eE/mctu„, B —> eB/mcw„. Here, 
X represents a vector in coordinate space, while cu„ is 
a chosen reference value for the normalising frequency, 
which, for instance, can be equal to the electron plasma 
frequency, or the frequency of the laser. In normalised 
units, the equations of motion of the particles are free 
of physical constants, except for a dimensionless coeffi¬ 
cient 


3mc^ 


( 10 ) 


which appears in the radiation reaction force term. 
Therefore, when including the radiation reaction force, 
the particle motion is no longer dependent only on the 
charge-to-mass ratio as in the Lorentz force E/,. This 
poses an additional challenge for the PIC implementa¬ 
tion of classical radiation reaction. 

In standard PIC codes, the system is represented us¬ 
ing macro-particles. Every macro-particle has the same 
charge-to-mass ratio as a single particle, and therefore 
the dynamics of the macro-particles is the same as the 
dynamics of the original particle species. However, af¬ 
ter including the radiation reaction, this no longer holds. 
If we examine equations (|^, Q and (0, we can see that 
the ratio between the radiation reaction force and the 
Lorentz force is proportional to the cube of the charge, 
and to the reciprocal value of the squared mass 


Frr ^ 

Fr nf- 


( 11 ) 


Let us consider a macro particle that represents tj 
electrons. The charge of the macro particle is = rje, 
and the mass is m„, = rim,,. Lor a single particle with the 
same mass and charge as the macro particle, the radia¬ 
tion reaction would be rj times stronger than in the case 
of a single electron: 


Frr {rjef 
Fl (rinieY mj 


( 12 ) 


and the trajectory of such particle would be different 
than the trajectory of a single electron (Lig. [^. This re¬ 
sult would be equivalent to assuming that rj electrons are 


radiating coherently. As a consequence, the results of a 
PIC simulation would be qualitatively different for dif¬ 
ferent number of particles per cell or different cell sizes. 
To obtain the correct dynamics of a macro-particle, it 
is therefore essential to use the real charge and mass 
to calculate the correct radiation reaction coefficient for 
a particular particle species. This approach yields the 
same result regardless of the macro-particle weight. 

3. Comparison of the models with standard radia¬ 
tion mechanisms 

To examine the physics captured by the different 
models, we compare the dynamics of a single particle 
with well known examples for which the particle trajec¬ 
tory and radiated power are known. 

We first consider synchrotron radiation, where a par¬ 
ticle moves in a constant external magnetic field. Taking 
only the Lorentz force into account, we expect the tra¬ 
jectory to be a perfect circumference due to the v x B 
term. When the particle has a very high initial mo¬ 
mentum, and is moving in an intense magnetic field, it 
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Figure 1: a) Trajectory of a macro-particle without radiation reaction 
b)-c) Trajectories of macro particles with weights t] = \ and t] = 5 
respectively, counter propagating with a laser pulse where ao = 100. 
Particle initial momentum is po = - 100. d) Energy of the same macro 
particles from a)-c) as a function of time. 
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LL 
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dt 
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ymc 


(9) 


Table 1: Radiation reaction contribution to the equations of motion: Fl - Lorentz force; p, e, m - particle momentum, charge and mass, y - 
relativistic factor; E, B - electromagnetic fields, c - speed of light, t - time. 


radiates, 
by HI: 


The radiative energy loss over time is given 


d^ 

dt 


2e^B^ 

3m‘^c^ 


(r 




(13) 


where ^ = ymc^ is the total particle energy and B is 
the magnetic field. If the particle is not relativistic, then 
^ ss mc^ and the right-hand side of ( |T3] l is close to zero, 
and therefore the energy loss is negligible. This is also 
the case when the B field is small. However, an ener¬ 
getic particle in a high intensity magnetic field will lose 
a significant amount of energy over time. As the parti¬ 
cle loses energy, its velocity decreases, and so does the 
curvature radius of its trajectory. This means that we do 
not expect a purely circular motion, but an inward spiral 
trajectory. When it loses a sufficient amount of energy, 
so that the right-hand side of ( [T3] l approaches zero i.e. 
^ ^ the particle trajectory then converges to a 

circumference. 

We have performed simulations for the same initial 
conditions using all the considered models (eqs. 0- 
(0). A typical example for a very strong damping is 


shown in Fig. [^where we can see that all the energy loss 
predictions closely resemble each other and match the 
analytical expression ( [T3| ). The strongest difference to 
eq. ( [T3| ) can be seen for the model BOS lH, which is still 
smaller than 0.3%. Therefore, we can confidently state 
that this effect is well resolved by all the models. This 
is not surprising because in this configuration, where B 
is constant and B ± p, equations 0-0 reduce to: 


BOS 

(±) ^-kB^t^p 

V dr Jrr y 


LL, S09, LLR, 
HOS, F93 

m ^-kB^p^p. 

V dr Jrr y ^ 

(14) 


For a highly relativistic particle, p » 1 and + \ ^ p^, 
so these two equations ([T^ are expected to yield similar 
results. 

Another illustration of the role of radiation reaction 
is Bremsstrahlung radiation. We let a charged particle 
propagate (yq = 10) along the x-axis starting from the 
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Figure 3: a) Electric field as felt by the pailicle over time; b) total 
pai'ticle energy vs. time; c) the difference in energy compared to the 
case without radiation reaction for k=0.001. 


Figure 2: Synchrotron radiation: electron trajectories and energy 
evolution. Tg - electron gyro period. The normalising frequency is 
the laser frequency corresponding to d = \fim. The initial conditions 
are: a) po = 100, Bo = 100, b) po = 100, Bq = 1000, c) po = 10, 
Bq = 5975. Dotted circle denotes the limit of the spiral motion. 
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origin in a static external electric field parallel to the 
particle motion given by (in normalised units): 

Eq I g(i-~-t)/a _ g-(L~x)la \ 

2 + e-iL~x)la j 

where Eq - a - A and L - 30 (see Fig. [^. The 
dominant term in eqs. 0-0 that leads to essentially 
the same result for the synchrotron radiation ( [l4| l, is 
identically equal to zero in this specific configuration. 
This allows us to explore the conceptual differences be¬ 
tween some of the models. In Fig. [^we observe that 
the LL model and the F93 model depend on the gradi¬ 
ent of the applied electric field. There is no radiation 
reaction observed in the BOS and the LLR models. The 
S09 and H08 models show some energy loss, but the 
rate at which the particle loses energy is constant where 
the electric field is constant and is not affected by the 
field gradient. 

This is better understood if we reduce the equations 
in Table 1 to the case where E||p: 

BOS and LLR (f= 0 
LLandF93 = + Ef) 

(16) 

[dtjRR- '^y+kpE 

H08 (^) 

\dt/RR ’ p- 

We see that the LL and F93 models reduce to the 
same expression. Their radiation reaction force depends 
only on the field derivatives, not on the absolute value 
of the field. There is no effect associated with the ra¬ 
diation reaction force in the BOS or in the LLR models. 
The S09 and HOS models lead to expressions that do not 
depend on the field derivatives, but only on the absolute 
value of the field, as can be observed in Fig. 

This analysis shows that if we choose one of the mod¬ 
els that are insensitive to the rate of the field change, we 
can apply it exclusively in a regime where this change is 
small enough not to significantly influence the motion. 
The condition that must be satisfied is then: 


(Fig. 1^) in a setup similar as in BtII . Here Eq stands for 
the Schwinger limit {Eq - irP'c^ f eh), which marks the 
transition from the classical to the quantum regime. For 
a field amplitude of Eq - 0.2 Ec the difference in the en¬ 
ergy loss predicted by different models is smaller than 
0.1% compared with the total electron energy (Fig. &)■ 
Since at E ^ Ec the classical models cannot be applied, 
a field amplitude outside the scope of applicability of 
the models is necessary to observe a non-negligible dif¬ 
ference between the models. 

Additionally, we would like to stress that the appli¬ 
cability condition ([T^ is satisfied in most physical sce¬ 
narios of interest. In particular, for laser pulses, the in¬ 
equality ( [T7] i is satisfied if aoy » 1, which is true when¬ 
ever radiation reaction is significant. In HSi this was 
confirmed by comparing the contribution of the particle 
spin and the contribution of radiation reaction force aris¬ 
ing from c/B/c/f and dEldt in the plane wave scenario. 
In this comparison, the spin gives a bigger contribution. 
In the classical regime that we are addressing here, the 
spin contribution is negligible, and so are the contribu¬ 
tions of dBIdt and dEldt. 

In summary, any of the proposed models can be used 
to describe the classical radiation reaction dominated 
regime. 

4. Testing the role of classical radiation reaction 
with the dynamics of electrons in intense laser 
pulses 

In order to examine the role of classical radiation re¬ 
action in scenarios with intense laser pulses and to deter¬ 
mine the conditions where such models should be used 
we consider the dynamics of a single electron interact¬ 
ing with a laser pulse. This is one of the main scenar¬ 
ios where radiation reaction can be explored ll50l ISTl 
and of very high relevance for future laser facilities 
||2l l45l l48l l52l l53l l54l l55]l . The laser pulse normalised 
vector potential is written as 15^ : 

A(x, t) - aQf{(p) cos0 e,, (IS) 


P 

trdc^ 


|FxP » 


dF 
dt ’ 


(17) 


where p is the particle momentum, is the perpendic¬ 
ular component of the Lorentz force F^ with respect to 

P 

To show that the contribution of the longitudinal 
fields to the radiation reaction in different models will be 
small in the vast majority of scenarios, we have chosen 
to illustrate a case with a very strong field Eq - 0.2 Ec 


where the temporal envelope f((f>) is a slowly varying 
function relative to the laser cycle {df I dt « cjq/), and 
the phase of the wave is given hy (f> - aiQt - wqx/c. We 
choose for /(0) a polynomial function IOt^ - 15t"^ + 6t^ 
where t = (pficoQ fpwHM) takes values in the domain 
[0,1] to define the envelope rise for 0 < i^/wo < tf^HM- 
Here, the pulse duration tfwhm is defined as full-width- 
at-half-maximum in the laser fields, and the envelope 
profile is symmetric relative to the point of the peak 
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intensity located at 0 = wq tpwHM- The relativisti- 
cally invariant normalised vector potential qq can be 
related to the linearly polarised laser intensity through 
ao = 0.8 W/cm2]/l[jum]. 

In the field of a linearly polarized laser pulse, a 
charged particle undergoes quiver motion. Without the 
radiation reaction force, the total particle energy re¬ 
mains unchanged after the interaction with the laser 
pulse. With radiation reaction, the total energy can de¬ 
crease substantially. If we consider a circularly polar¬ 
ized laser pulse, the situation is similar. In Fig. |^this is 
illustrated for a linearly and a circularly polarized laser 
with identical temporal envelopes and the same intensi¬ 
ties = 100, = 100/ V2, do = Ijum). The ini¬ 

tial normalised momentum of the particle is po - 100, 
opposite to the laser propagation direction. The total 
energy that the particle loses while interacting with the 
laser is about 30% and is the same for the linearly and 
the circularly polarized case. All the models give sim¬ 
ilar results (the biggest difference in the final energy in 
Fig.gis 0.03%). 

It is important to estimate the threshold when radia¬ 
tion reaction starts to significantly affect the motion. In 
Fig. B we see that for counter-propagating particle and 
a laser pulse there is already strong radiation reaction 
for ao - 100 if po - 100. However, if we repeat the 
simulation with the same laser parameters for a less en¬ 
ergetic particle, radiation reaction has a weaker effect. 
For a particle starting at rest, there is no radiation reac¬ 
tion observed for aq = 100. This indicates that in order 
to estimate whether radiation reaction is important for 
a set of particular conditions, it is also required to have 
the information about the relative motion of the particle 
and the laser. We can assume that the radiation reac¬ 



tion force is significant when the following condition is 
satisfied: 

^ > lO-l (19) 

lI 

In eq. 0, the dominant term due to radiation for high 
y (ultra-relativisic particle) written in units where all 
quantities are normalised to the laser frequency is: 


2 e^ojQ 

3 mc^ 


rp 



(E ■ p)" 


( 20 ) 


In these units E,B ao and y ~ p, and thus 
the left hand side of eq. is on the order of 

(2e^a)ol(3mc^))y^aQ for a particle counter-propagating 
with the wave. For a laser with a wavelength do ~ Ipm 
eq. ([T9]l then becomes: 



The limits of validity of the classical approach for ra¬ 
diation reaction have been discussed in 15311521. When 
the electric field of the laser approaches the Schwinger 
field in the rest frame of the particle, quantum effects are 
expected to dominate. For a head-on collision, where 
radiation reaction has the strongest effect on the motion, 
the classical equations can then be used if 



Therefore, from eq. ( |2T] i and eq. 122) 1 we can identify 
the range of applicability of the classical radiation reac¬ 
tion models described here for the case of an electron 
colliding head-on with a laser as: 



For instance, a 300 J laser pulse with A - 1/im and 
30 fs duration, yields a peak power of 10 PW and, 
when focused to a 10 pm focal spot, a peak intensity 
of / = 10^^ W/cm^. The corresponding vector poten¬ 
tial is ao ~ 100. Future laser facilities will be able to 
provide normalised vector potentials of this magnitude, 
and for yo ^ 100 the inequalities (|2^ are verified. 


5. The role of the timestep for simulating particle 
motion in an intense electromagnetic wave 


Figure 4: Particle energy over time in the field of linearly (solid line) 
and circularly (dashed line) polarized laser pulse. 
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The dynamics of electrons in intense electromagnetic 
waves can be strongly relativistic. At ultra high in¬ 
tensities the relative dynamics of the particles in the 







































laser field is highly nonlinear and the relevant numer¬ 
ical parameters must be selected with extra care. In this 
section, we will see that the choice of the simulation 
timestep is very important for the accuracy of the final 
result, and that the optimal value depends on the laser 
intensity, electron energy, total simulated time, as well 
as whether radiation reaction is significant or not. 

The higher resolution requirements at ao » 1 were 
reported previously in fSSl . where they studied the par¬ 
ticle motion in an intense electromagnetic plane wave in 
ID geometry and found the approximate convergence 
condition cAf/d « I/aq (here A and aq are the laser 
wavelength and the normalised vector potential, c is the 
speed of light and Af the simulation time step). Their 
study focuses on long interaction at relativistic intensi¬ 
ties with Ao ~ 10. The higher resolution requirement 
is even more relevant in scenarios where radiation re¬ 
action plays an important role because of the prospects 
of dealing with ultra-high intensities (aq » 10). How¬ 
ever, the state-of-the-art and near-future laser facilities 
that aim at reaching such intensities are usually short 
(~ 30 fs), so we opt to consider short laser duration in 
our further analysis. In fact, as discussed below, we will 
have different resolution requirements depending on the 
particle energy, momentum direction and the maximum 
intensity of the laser field. For our case, the resolution 
requirements could be less strict than that of ref. Il58l 
depending on the energy of the interacting particles. As 
the average energy loss due to radiation reaction also de¬ 
pends on the laser duration ll5^ . we could expect differ¬ 
ent resolution requirements for different laser durations 
in cases where there is significant radiation reaction. 

An initially small integration error can, in principle, 
grow during a long simulation. To assess this effect, we 
perform temporal resolution tests for realistic simula¬ 
tion timescales using the fourth order Runge-Kutta in¬ 
tegration method. Usually, the longest simulations that 
can be performed using PIC codes are the ones that can 
take advantage of the moving window technique used 
to simulate only a portion of plasma around the laser 
pulse 15^ . We therefore take a typical laser wakefield 
acceleration (LWFA) setup for a 0.5 GeV electron ac¬ 
celeration stage operating in matched conditions Il60l as 
an example of one such long simulation. For a driver 
laser wavelength Ald = 1 /tm, it is required to have 
a 1.6 mm plasma length, laser aqzj - 1-6 and the ra¬ 
tio between the laser and the electron plasma frequency 
oj() jupe - 10. This amounts to a total simulation time 
of Ttotai - 10"^ where the laser frequency is given 
by aiQ = 1.88 X lO'^ s“'. We take Ttotai to be a char¬ 
acteristic time of a long simulation. The conclusions 
drawn by using Ttotai will be applicable for any simula¬ 


tions where simulation time Tsim ^ TtotaU which applies 
to most PIC simulations performed with currently avail¬ 
able resources. 

Now that we established a typical simulation time, 
we consider the motion of an electron in a transversely 
plane electromagnetic wave with a longitudinal enve¬ 
lope that propagates along the x-direction. The elec¬ 
tromagnetic wave is expressed analytically according to 
the eq. ( [T8] l and has 30 fs duration at FWHM in inten¬ 
sity. The electron is considered as a test-particle that 
gives no feedback to the fields, and its motion is fol¬ 
lowed always for the same amount of time Ttotai- The 
electron initial momentum is chosen to be parallel to 
the wave propagation direction and either co- (po > 0) 
or counter-propagating (po < 0) with the laser. We have 
studied the convergence of the particle trajectories for 
different initial momenta po, and different vector poten¬ 
tial of the wave aq varying the timestep, while keeping 
the simulation duration in the units normalised to the 
laser frequency ojq equal to Ttotai. 

For some examples when po > 0 the total interaction 
time Tint of the particle with a laser may be much longer 
than the time we are interested in Ttotai- Hence, for our 
purpose the resolution required for a reliable simula¬ 
tion ensures correct numerical integration of trajectories 
of such particles within the simulation time Ttotai- We 
would like to warn the reader that this resolution may 
be insufficient if one wishes to simulate the entire dura¬ 
tion of laser-electron interaction if Tim » Ttotai ■ 

In order to establish a quantitative convergence crite¬ 
ria that would allow an automatic evaluation, we have 
first performed the simulations with a very small time 
step dt - 10“"^ well beyond the convergence limit 
for each individual case. Then, we chose 1000 refer¬ 
ence points equally spaced in time (for T=10, 20, ...) 
along the trajectory which would serve as a benchmark 
to compare with the runs with longer dt. We can then 
define the relative error as: 


Zr ■ 

Here, y is the Lorentz factor of the converged result, |Ay| 
stands for the absolute difference between this value and 
the one from the current test-particle run and the sums 
are over all reference points. We take the error to be ac¬ 
ceptable if R <1% for which we consider that the result 
has converged. An illustration of converged and uncon¬ 
verted particle trajectories for aq = 40 and po - -3 is 
shown in Fig. The figure shows that the condition 
R <1% is conservative because even a small visible de¬ 
viation in a particle trajectory leads to R »1%. 

Comprehensive tests were done for a wide range of 
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Figure 5: Illustration of converged and unconverted particle trajectories, a) Longitudinal coordinate over time, b) Transverse momentum over time, 
c) Longitudinal as a function of transverse momentum, d) Transverse coordinate over time, c) - d) Upper plots show the difference of the particle 
trajectory obtained with dt = 0.2 (dashed red line) compared with the converged case (solid black line). Middle plots show the same comparison 
for dt = 0.\ (blue dashed line), and lower plots for dt = 0.05 (green dashed line). 


values of qq and initial momenta. A subset of these is 
plotted in Fig. where each pixel represents a single 
simulation (300x100 pixels in each panel). It can be 
seen that below aq = 5, the condition of having ~ 30 
points per laser period is sufficient (here equivalent to 
dt-0.2 Wq'). If flo is higher, the most demanding case is 
when the particle has very low energy before the inter¬ 
action with the laser. In that case, the particle can get a 
strong kick, on the order of its total energy, within a sin¬ 
gle time step. This can modify its direction of motion 
and put the particle on a different trajectory, allowing 
the error to grow as the particle gains energy in the laser 
field. If a particle is already ultra-relativistic, the same 
error in the accelerating force does not change the mo¬ 


tion significantly. 

Our findings are consistent with the ones reported in 
Ref. Il58l . They have shown that the ’’stopping point” 
of electron trajectories is the most sensitive section with 
respect to integration errors. Resolution requirements 
are therefore expected to be more strict for particles that 
do have a point of zero-momentum during the laser in¬ 
teraction. 

From Fig. we can conclude that for studying the 
interaction of relativistic electron beams with intense 
lasers where y > ao/3, a resolution of 30 timesteps per 
laser period can still be kept. For smaller values of ini¬ 
tial 7 this no longer holds and special care should be 
applied when simulating intense lasers interacting with 
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Figure 6: Convergence map for different time steps at a) high and b) low intensities. Convergence map for c), d) dt = 0.05, q) dt = 0.1, f) 
dt = 0.025, with and without accounting for radiation reaction. Color legend for both a) and b) is presented in panel b), while for c) - f) it is 
presented in panel c). The laser pulse duration was 30 fs (==; 56.5 and total simulation time for all simulations Ttotai = iO^ Below 

ao = 100, the convergence does not depend on whether radiation reaction is included or not. Below aq = 5, the standard condition (dt=0.2 ~ 30 
points per laser period) is enough. For higher values, the most limiting case is laser in a cold plasma - the interaction of ultra relativistic particles 
can still be modelled with dt = 0.2. The difference in the regions with higher values of gq with and without radiation reaction can be attributed to 
the particle energy loss that lowers the effective p of the particle. 
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Figure 7: Particle trajectories for uq = 300 and po = -40. Transverse momentum p 2 versus time a) without and b) with radiation reaction. 
Thick black lines correspond to the fully converged trajectories. The green lines corresponds to dt = 0.05, for which the simulation has converged 
without radiation reaction and has not converged when radiation reaction is taken into account. This is even more evident when we compare all the 
trajectories on a smaller time interval near the end of the interaction c) without and d) with radiation reaction. 


cold plasmas. These are general conclusions that ap¬ 
ply whether or not the radiation reaction is accounted 
for in the code. We have repeated the same procedure 
with and without radiation reaction, and for lasers where 
flo < 80 we have not observed any differences in the 
convergence map (Fig. @a,b). 

A difference between the convergence map with and 
without RR (Fig. c-f) arises only for very high aq- 
This is not surprising since in these cases the particles 
lose significant energy due to the strong radiation re¬ 
action, and thus move to the regions in the parameter 
space of Fig. j^where smaller dt would be required. For 
instance, if at the start of the interaction with a laser of 
aq = 300, an electron has po = “40 me, after the inter¬ 
action it will have p = -15 me. Figure]^ shows that 
dt — 0.1 is not small enough for the particle trajectory 
to converge, with or without radiation reaction (the cor¬ 


responding trajectories are illustrated in Fig. [^coloured 
blue). The electron trajectory satisfies our convergence 
criteria without radiation reaction for dt - 0.05, but 
not if we take radiation reaction into account (Fig. igi 
and Fig. green). However, for po = -15, the con¬ 
vergence criteria without radiation reaction is also not 
satished for dt=0.05 because lower particle momentum 
requires higher resolution (it is easier to reach the ’’stop¬ 
ping point” with a lower momentum). This means that 
by causing energy loss, and thus lowering the particle 
momentum, radiation reaction leads to more stringent 
convergence criteria that is evident in Fig. |^c-f. If we 
consider the lowest value of the average momentum in 
the simulation (i.e. the final momentum), we can apply 
the convergence map obtained without radiation reac¬ 
tion. A way to quickly estimate the total energy loss dur¬ 
ing a head-on collision with a laser pulse can be found 
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- = cV X B — 47rj 


f = -cVxE 


Model 

Overhead 

Added quantities 

BOS 

1 

0 

LL 

4.5 

6 per grid point 

S09 

2 

0 

HOS 

1.5 

3 per particle 

F93 

2 

3 per particle 

LLR 

2 

0 


Table 2: Computational cost for each RR model 


Figure 8: Main PIC loop, x,, u/, F/ - particle position, generalized 
velocity u/ = yv, and force acting on the i-th particle; j-total current; 
E, B the electric and the magnetic field respectively. 

in Ref. l52l- 

6. Computational overhead of radiation reaction 

modelling 

To assess the computational overhead of the radiation 
reaction schemes for particle-in-cell codes, we need to 
consider the full structure of the PIC loop (presented 
in Fig. 1^ and the main aspects that affect its perfor¬ 
mance. The electromagnetic fields in the simulation are 
stored on the grid nodes, while particles explore the full 
6D phasespace. The helds stored on the grid are inter¬ 
polated to the particle positions and then the particles 
are pushed according to the Lorentz force. Later, the 
currents created by all the particles are deposited on the 
grid, where this information allows to advance the fields 
on the grid nodes. The new fields are then interpolated 
again to the new particle positions to start the next it¬ 
eration. OSIRIS is parallelised through a spatial do¬ 
main decomposition, where only the neighbouring com¬ 
puter nodes are required to exchange information. This 
exchange consist of sending the information about the 
particles that moved to the sector of space that belongs 
to the neighbouring node and about the electromagnetic 
fields in the vicinity of the boundary between the two. 

Introducing the radiation reaction into the PIC loop, 
in principle, does not affect the parallelisation and 
should not affect the overall scalability of the PIC code. 
The main difference is in the particle pusher, where in¬ 
stead of the Lorenz force, it is necessary to apply a full 
equation of motion that also incorporates the radiation 
reaction force. The particle pusher in OSIRIS can fol¬ 
low a second-order Boris method ED, or a fourth-order 


Runge-Kutta integration method for the particle equa¬ 
tions of motion. Radiation reaction can be implemented 
in both of them (see ref. l48l for instructions on adapt¬ 
ing the LL method for Boris pusher). 

When examining the equations Q-(|^ one can expect 
a significant overhead of the particle pusher as com¬ 
pared to the case where we only had the Lorentz force. 
To estimate the additional computational cost, we have 
determined the particle pusher overhead for every model 
used here. The overhead is dehned as the ratio between 
the number of floating point operations (flop) required 
for the radiation reaction force and for the Lorentz force 
respectively. 

In state-of-the-art supercomputers, memory access is 
often more time-consuming than additional flop. The 
balance depends on a particular computer architecture, 
but in all cases it is advisable to minimize the amount of 
saved quantities. The computation of radiation reaction 
force requires sometimes additional memory. For ex¬ 
ample, in eq. the Lorentz force from a previous time 
step needs to be saved. This implies that three space 
components of the force are to be added for each particle 
in the simulation. Also, the LL model requires to save 
electric and magnetic helds from the previous timestep, 
which results in six extra quantities per grid point. In 
addition to adding memory requirements (more RAM 
needed for the problem of the same size), this would 
also double the amount of communications related to 
the exchange of held values near the boundary between 
the neighbouring nodes. 

The additional computational cost is summarized in 
Table 2. The models that do not require any additional 
memory (BOS, S09 and LLR) also have low overhead in 
terms of Hop. 

Therefore, any of these three models can be intro¬ 
duced without affecting signihcantly the performance of 
a PIC code. We chose the LLR model and implemented 
it in OSIRIS pusher with the fourth-order Runge Kutta 
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integration method. 

7. Conclusions 

We have compared different approaches to include 
classical radiation reaction in a PIC code. Out of the 
numerous models available in the literature, we chose a 
subset compatible with the standard PIC algorithm and 
compared their performance in modelling the electron 
motion in physical scenarios of interest. 

The first test was on the electron motion in a con¬ 
stant B-field (i.e. synchrotron), where all the models 
gave the same results, with a minimal departure of BOS 
at extremely high fields {B ^ 0.1 Be). Interaction of 
electrons with a counter-propagating laser pulse (both 
for circular and linear polarisation) also did not show 
any appreciable distinction between the models. This 
is not surprising, because for relativistic electrons the 
leading order term in the radiation reaction force asso¬ 
ciated with the transverse acceleration is the same for all 
the compared models. However, differences arise in the 
presence of a longitudinal electric field with a spatio- 
temporal gradient. For such fields, the model F93 pre¬ 
dicts the same recoil as the LL model, that depends on 
the electric field gradient. The models S09 and H08 de¬ 
part from the LL results, but do predict an energy loss 
that depends solely on the field magnitude (not the gra¬ 
dient). The models BOS and LLR do not account for 
any energy loss in this configuration. Even though this 
offers a qualitative distinction between the models, the 
electric field needs to be on the order of the Schwinger 
critical field to show it, and even then the differences 
make for less than a percent of the total electron energy. 
They are, therefore, unlikely to make any real impact to 
the simulation results for realistic configurations. For 
example, these differences are negligible compared to 
the leading order term of the radiation reaction force for 
particles in the field of a laser pulse if » 1, which 
is necessarily satisfied in case of a significant radiation 
reaction. 

One additional question that needed to be addressed 
for PIC implementation is the definition of radiation re¬ 
action for a macro-particle. We note that the radiation 
reaction depends nonlinearly on charge and mass. To 
avoid overestimating the output radiation by artificially 
assuming coherent emission, the real charge and mass 
of a single particle should be used when calculating the 
radiation reaction force for a macro-particle. 

The models are all fit to describe the physics, but the 
computational cost is minimal (in terms of flop and ad¬ 
ditional memory required) for the models BOS, S09 and 
LLR. They can be included in a massively parallel PIC 


code just by re-writing the particle pusher, without any 
additional quantities to store. We have opted to imple¬ 
ment the LLR model in OSIRIS. 

Apart from choosing a model that accounts for radi¬ 
ation reaction, to correctly model the particle dynamics 
at extreme intensities it is essential to use the appropri¬ 
ate temporal resolution in the simulations. For the same 
laser frequency, a higher ao requires higher resolution to 
account correctly for the electron motion (especially if 
the electron is not relativistic at the start of the interac¬ 
tion). Detailed analysis supported by test-particle simu¬ 
lations was performed to facilitate the choice of tempo¬ 
ral resolution for simulations at high intensities. 

Taking all the aforementioned into account, it is pos¬ 
sible to simulate the classical radiation reaction domi¬ 
nated regime with PIC codes by keeping the same paral¬ 
lel structure and introducing modifications only at parti¬ 
cle push-time. Such modifications should not affect sig¬ 
nificantly the computational performance of massively 
parallel PIC simulations, with several models (e. g. 
LLR or BOS) capturing the appropriate dynamics with 
minimal computational overhead. 
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